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Certainly  the  most  useful  background  reference  for  finite  element 
modelling  of  reinforced  concrete  is  the  ASCE  Task  Committee  Report  [6]. 
This  publication  represents  the  state-of-the-art  through  1982. 

A  recent  review  article  on  the  finite  element  analysis  of  reinforced 
concrete  shell  structures  describes  the  current  state-of-the-art  (Ramm  [19]).  A 
main  issue  seems  to  be  the  competition  between  classical  elastic  buckling 
effects  and  strength,  or  material  modelling,  effects.  Most  reinforced  concrete 
shells  in  service  are  very  thin  compared  with  structures  of  interest  to  NCEL. 
Thus  some  issues  of  interest  to  the  reinforced-concrete  shell  community  are 
not  germane.  The  relatively  thick  slab-like  magazine  structures,  with 
important  three-dimensional  built-up  sections  (e.g.,  around  columns),  of 
interest  to  NCEL  make  elastic  buckling  theory  largely  irrelevant.  Issues  of 
material  modelling  thus  become  of  prime  interest  as  do  large-deformation 
factors  due  to  the  intensity  of  design  loadings  and  concern  with  ultimate 
strength.  Most  currently  employed  procedures  use  typical  assumptions  in  the 
development  of  shell  elements  and  there  seems  to  be  a  preference  for 
"degenerated,"  or  continuum-based,  shell  elements.  Plane  stress  constitutive 
algorithms  are  emphasized  (see,  e.g.,  Ramm  [19,  20])  and  variants  of  the 
Kupfer  bi-axial  failure  envelope  [16]  are  typically  used.  In  compression  and 
tension,  failure  mechanisms  of  crushing  and  cracking,  respectively,  are 
assumed.  Both  involve  strain-softening  which  is  perhaps  the  most 
important,  delicate,  and  controversial  aspect  of  concrete  modelling.  Softening 
due  to  cracking  is  often  referred  to  as  tension  stiffening.  Some  form  of 
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smearing  is  usually  employed  to  model  cracking.  Typically,  elastic  moduli  are 
reduced,  or  set  to  zero,  to  manifest  the  presence  of  a  one  or  more  cracks  at  a 
point.  Cracks  may  be  assumed  to  have  fixed  orientation,  or  rotate  with 
principal  stresses  [1,  13].  This  is  also  a  controversial  aspect  of  concrete 
modelling  raising  issues  of  invariance  and  micro-mechanical  physical 
mechanisms  [2-4].  To  account  for  rough  crack  surfaces  and  doweling  some 
form  of  shear  stiffness  retention  factor  is  frequently  employed.  Tension 
stiffening  is  often  assumed  to  subsume  phenomena  associated  with  concrete- 
reinforcement  bond  slip.  This  is  obviously  rather  crude  at  best. 

Recent  work  of  Cervera  et  al.  [7,8]  is  rather  advanced  in  many  respects 
and  appears  particularly  germane  to  the  interests  of  NCEL.  Cervera  et  al.  are 
concerned  with  three-dimensional  states  of  stress  and  especially  transverse 
shear  states  in  the  vicinity  of  supports.  They  eschew  shell  elements  in  favor 
of  three-dimensional  finite  element  models,  in  particular,  8-20  node 
isoparametric  bricks.  They  argue  against  reduced  integration  procedures  and 
reassess  the  range  of  applicability  of  brick  elements  in  bending  situations. 
This  imposes  a  length-thickness  ratio  limit  in  modelling  which  may  be 
extreme  depending  on  the  order  of  element  employed.  They  model 
reinforcement  by  way  of  perfectly  bonded  membranes  having  unidirectional 
properties  "equivalent"  to  the  distribution  of  steel  reinforcing  bars.  The 
membrane  and  concrete  brick  elements  are  assembled  together  in  the  usual 
fashion.  The  assumption  of  perfect  bond  is  a  shortcoming,  nevertheless  it  is 
typical  of  the  state-of-the-art  in  large-scale  structural  modelling.  The  failure 
envelope  in  compression  is  given  by 

/  (I,,  I2)=  <«l,  +  3PI2)'12  =  a0 
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where  /  2  and  /  2  are  the  first  and  second  stress  invariants,  respectively, 
and  a  and  /3  are  parameters.  These  can  be  used  to  fit  Kupfer's  bi-axial 
data.  An  associative  flow  rule  is  assumed.  Hardening  is  assumed  to  decrease 
in  this  model  eventually  becoming  a  perfectly-plastic  model.  A  strain-defined 
crushing  condition,  when  satisfied,  indicates  release  of  all  stresses  and 
stiffness.  This,  of  course,  is  a  form  of  strain  softening.  Two  fixed  cracks  are 
allowed  to  form  at  each  sampling  point.  After  the  first  crack  is  formed,  the 
second  crack  may  only  form  in  the  plane  perpendicular  to  the  first  and  the 
formation  of  the  second  crack  depends  only  on  the  two-dimensional  state  of 
stress  perpendicular  to  the  plane  of  the  first.  A  cracking  failure  surface  in 
stress  space  represents  a  three-dimensional  generalization  of  the  Kupfer-type 
bi-axial  surface  in  tension.  Strain-softening  in  tension  is  assumed  to  obey 
relations  emanating  from  fracture  mechanics  which  engender  mesh- 
dependent  softening  moduli  which  in  turn  desensitize  results  to  the  degree  of 
mesh  refinement.  These  ideas  are  described  in  more  detail  later.  It  is 
mentioned  that  tension  stiffening  due  to  reinforcing  may  be  "heuristically 
included  by  assuming  a  higher  fracture  energy  (release  rate)  for  reinforced 
concrete  than  for  plain  concrete."  A  variable  shear  retention  factor  is 
employed.  Cracking  is  assumed  to  reduce  the  compressive  hardening 
modulus  in  transverse  directions.  Static  as  well  as  dynamic  studies  have  been 
performed  by  Cervera  et  al.  [7,  8]. 


Cracking  and  crushing  of  concrete  are  frequently  modelled  with  strain- 
softening  mechanisms.  From  a  phenomenological  point  of  view,  there  is  no 
doubt  that  softening  occurs.  However,  it  is  often  and  convincingly  argued 
that  softening  is  not  a  constitutive  phenomenon,  but  rather  a  structural 
phenomenon  indicative  of  progressive  damage.  Consider  a  concrete  tensile 
specimen.  In  a  displacement-controlled  test,  it  is  observed  experimentally 
that  a  peak  tensile  force  is  attained  and  subsequently  the  force  drops  off 
toward  zero.  This  softening  is  accompanied  by  a  localized  band  of 
microcracking  and  crack  coalescence.  In  the  undamaged  region  of  the 
specimen,  elastic  unloading  takes  place.  If  the  softening  branch  of  the  force- 
displacement  diagram  is  interpreted  as  a  material  property  and  used  to 
identify  a  softening  modulus  in  a  constitutive  equation,  pathological  results 
can  ensue.  To  see  this  in  as  simple  a  setting  as  possible,  assume  the  softening 
branch  is  linear  and  the  one-dimensional  stress-strain  response  is  as  indicated 
in  Figure  1.  The  notation  is  as  follows: 

cr  =  Tensile  stress 

oy  =  Value  of  c  at  which  fracturing  commences 

£  =  Tensile  strain 

£j  =  Value  of  e  at  which  fracturing  commences 

£q  =  Value  of  £  at  which  specimen  ceases  to  carry  any  load 

£  =  Young's  modulus 

£s  =  Softening  modulus 
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It  is  perhaps  apparent  that  the  strain  computed  by  dividing  the  specimen 
displacement  by  the  specimen  length  in  no  way  represents  the  strain  of  any 
material  point  in  the  specimen.  The  major  portion  of  the  specimen  is 
unstrained  due  to  elastic  unloading  whereas  a  small  portion  has  fractured, 
which  may  be  interpreted  as  infinite  strain.  Thus  employing  an  average 
strain  in  the  development  of  a  constitutive  equation  clearly  is  problematic. 
Nevertheless,  this  process  has  been  frequently  employed  in  practice.  Consider 
a  one-dimensional  finite  element  model  of  the  experiment.  Assume  a 
discretization  of  piecewise  linear  elements  of  length  h  .  The  constitutive 
behavior  of  each  element  is  governed  by  Figure  1.  No  matter  how  many 
elements  are  employed,  the  stable  solution  to  monotonically  increasing 
applied  displacement  is  for  one  element  to  follow  the  stress-strain  diagram  of 
Figure  1  and  all  others  to  elastically  unload  as  soon  as  £  exceeds  ty.  This 
can  be  seen  by  imagining  that  the  values  of  oy  for  each  element  are 
randomly  perturbed.  In  this  case  one  element  must  soften  before  the  others 

and  by  virtue  of  equilibrium  the  others  must  experience  a  decrease  in  stress 
prior  to  £  exceeding  their  fy's.  Clearly  the  solution  to  the  unperturbed 

problem  is  non-unique.  How  well  does  the  solution  represent  the 
experiment  which  was  used  to  derive  the  constitutive  behavior?  Let  us 
compute  the  energy  dissipated  in  the  calculation.  This  is  given  by  the  area 
under  the  stress-strain  diagram  multiplied  by  the  volume  of  the  element.  Let 
A  =  Cross-sectional  area  of  undamaged  specimen. 

Then  the  energy  dissipated  is 
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Note  that  this  expression  is  O  (h  ).  Consequently,  the  results  exhibit  a  mesh 
sensitivity  which  is  totally  spurious.  The  only  way  the  results  of  the 
experiment  can  be  reproduced  is  if  h  is  taken  equal  to  the  specimen  length. 

This  obvious  pathology  needs  to  be  circumvented  otherwise  numerical 
calculations  become  meaningless.  A  number  of  authors  have  proposed 
viewing  the  softening  modulus,  £s  ,  as  a  function  of  the  mesh  so  that 

numerical  calculations  become  desensitized.  Representative  of  these  works 
are  studies  of  Wiliam  and  his  colleagues  (see  [22  -  24]).  The  softening 
modulus  is  determined  so  that  the  energy  dissipated,  (1),  remains  constant, 
viz., 


h  A  a 
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(l  1  } 

{E+Es(h)  J 


-  A  ■  Gj  -  constant 


(2) 


where 

Gj  =  Fracture  energy  release  rate. 

From  (2), 

Es(h)=Y—  (3) 

2 

where  L  =  2  G  /  Gf  denotes  a  characteristic  scale  associated  with  the 


damaged  region. 


Remarks 
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2.  Wiliam  et  al.  [22]  assert  that  (3)  imposes  a  maximum  element  length, 
namely  h  <,  C.  Although  there  seems  to  be  no  mathematical 
problem  with  h  >  t  (see  Figure  2),  typical  "strain  driven"  constitutive 
algorithms  would  not  be  able  to  handle  this  situation. 

3.  Frequently,  tensile  cracking  is  modelled  with  strain-softening  elastic- 
plastic  constitutive  equations.  It  is  interesting  to  calculate  the  plastic 
modulus,  H,  corresponding  to  the  softening  modulus  given  by  (3).  By 


definition 


and  so 


-EH 

Es  ~  E  +  H 


H  ( h)=~~E 


The  mesh-dependence  of  H  is  apparent.  It  appears  that  this  simple 
way  of  addressing  the  mesh  sensitivity  problem  was  first  proposed  by 
Pietruszczak  and  Mroz  [18]  who  derived  (5)  from  an  entirely  different 
point  of  view.  See  also  Zimmermann  [25]  and  Wiliam  et  al.  [23]. 
Similar  ideas  have  been  introduced  by  Bazant  (see,  e.g.,  Bazant  and 
Cedolin  [5]),  Hillerborg  et  al.  [14],  and  Nilsson  and  Oldenburg  [17]. 

4.  The  generalization  of  (5)  to  multidimensional  applications  is  done  in 
ad  hoc  manner.  Various  authors  have  employed  (5)  with  h  being 
taken  as  simply  the  square-root  of  the  area  of  the  element  in  two 
dimensions  and,  likewise,  the  cube-root  of  the  volume  in  three 
dimensions  (see,  e.g.,  Zimmermann  [25]).  For  this  strategy  to  be 
reasonable  elements  need  be  as  close  to  squares  or  cubes  as  possible. 
Elements  with  large  aspect  ratios  clearly  necessitate  h  being  defined 


f  . ■»%.« \1 
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adaptively  as  a  characteristic  element  dimension  perpendicular  to  the 
direction  of  maximum  principal  stress.  When  higher-order,  multiple- 
quadrature-point  elements  are  employed  a  subelement  area  or  volume 
may  be  assigned  to  each  quadrature  point  for  calculating  local  values  of 
h  (see,  e.g.,  Cervera  et  al.  [7, 8]). 

5.  Modelling  fracture  damage  within  the  formal  structure  of  classical 
plasticity  theories  does  not  seem  consistent  with  the  observation  that 
unloading  paths  do  not  in  general  coincide  with  those  governed  by  the 
elastic  moduli.  Introduction  of  elastic  damage  mechanisms  causes 
unloading  at  reduced  values  of  the  elastic  modulus.  A  pure  elastic 
damage  model  unloads  towards  the  origin  (see  Fig.  3).  The  presence  of 
softening  creates  the  same  potential  for  mesh  sensitivity  within  the 
context  of  elastic  damage  mechanisms  as  described  previously.  This 
should  be  obvious  because  the  illustrated  pathology  has  nothing 
whatsoever  to  do  with  unloading.  Recently  Resende  [21]  has  proposed 
an  isotropic  damage  theory  for  concrete.  Despite  the  presence  of 
softening,  no  attempt  is  made  in  this  formulation  to  desensitize  the 
mesh  dependency.  Numerical  results  presented  in  [21]  clearly  exhibit 
this  phenomenon.  It  is  of  course  possible  to  desensitize  an  elastic 
damage  model  in  the  same  way  as  described  previously.  For  example, 
consider  the  following  one-dimensional  model: 

G  =  (l  -d)Ee  (6) 

where  d€  [0,1]  is  the  damage  parameter.  The  accumulation  of  damage 
is  governed  by  the  damage  evolution  law: 
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Where  x  =  ^E£z  and  h  is  the  damage  evolution  function.  (A 
more  detailed  description  of  theories  such  as  this  may  be  found  in  [15].) 


•  /j  . 

d  =  -  Eee 

X 


(7) 


Taking  the  time  derivative  of  (6)  and  employing  (7)  yields 


o  = 


(l-d)E--(Ee)2 

x 


(8) 


The  expression  in  brackets  is  the  tangent  modulus  during  damage 
loading.  At  this  point  h  is  not  specified.  It  may  be  determined  so  that 
the  softening  branch  of  the  stress-strain  diagram  is  reproduced  as 
before,  viz., 

A  -F. 

(9) 
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When  a  combination  of  damage  mechanisms  is  employed  in  a  theory, 
such  as  in  an  elastoplastic  damage  theory  [15],  it  is  not  so  clear  how  to 
perform  the  necessary  modifications  to  desensitize  mesh  dependency. 
The  ad  hoc  procedure  of  introducing  the  mesh  parameter  h  into 
constitutive  equations  thus  does  not  seem  satisfactory  in  more  complex 
situations. 

6.  Softening  branches  of  stress-strain  diagrams  are  frequently  modelled 
with  more  complex  functions  (e.g.,  exponentials,  see  Cervera  et  al.  [7,  8], 
Nilsson  and  Oldenburg  [17]  ).  This  poses  no  essential  problems. 


Due  to  turning-point  phenomena  in  solutions  of  quasi-static 
reinforced-concrete  problems,  it  is  necessary  to  employ  "arc-length"  strategies. 
(These  are  referred  to  as  "continuation  methods"  in  the  mathematics 
literature.)  Crisfield  has  been  a  leader  in  the  development  and  application  of 
these  methods  to  reinforced-concrete  structures.  Representative  papers  are  [9- 
12].  Any  computer  program  for  the  analysis  of  reinforced-concrete  structures 
should  be  equipped  with  algorithms  of  this  type.  Unfortunately,  even  the 
best  procedures  encounter  difficulties  when  faced  with  severe  bifurcation  and 
unloading  phenomena  primarily  brought  about  by  softening  mechanisms 
present  in  reinforced-concrete  constitutive  models. 

Most  applications  of  current  interest  to  NCEL  involve  dynamic 
phenomena  (e.g.,  magazines  subjected  to  blast  wave  loading).  Thus  state-of- 
the-art  implicit  and  explicit  techniques  are  required.  In  addition,  due  to  the 
complexity  of  reinforced-concrete  behavior,  automatic  control  of  time-step 
size  and  nonlinear  solution  strategy  within  each  step,  etc.,  should  be  included 
in  a  computer  program.  It  is  felt  that  a  much  greater  degree  of  algorithmic 
automation  than  present  in  existing  techniques  (see,  e.g.,  Cervera  et  al.  [7,  8]  ) 
will  be  desired  and,  in  fact,  required  for  the  reliable  and  efficient  solution  of 
problems  anticipated. 


Steel  reinforcement  is  modelled  as  an  elastoplastic  material  in 
traditional  fashion.  The  most  vexing  issue  in  reinforcement  modelling  is  the 
bond  law  governing  the  steel-concrete  interface.  Studies  have  been 
performed  in  which  discrete,  nodally  assembled  reinforcing  bars  have  been 
modelled  and  bond-slip  laws  assumed.  These  studies  have  typically  been  of 
simple  experimental  configurations  in  which,  at  most,  only  a  few  bars  are 
present.  Chapter  3  of  [6]  contains  a  good  summary  of  work  through  1982. 
Complex  reinforcing  patterns  in  engineering  structures  have  been  modelled 
by  either  distributed  and/or  embedded  equivalent  reinforcing  elements. 
Perfect  bond  is  typically  assumed.  This  important  aspect  of  the  overall  failure 
strength  thus  is  neglected  by  state-of-the-art  approaches  and  thus  must  be 
considered  a  serious  deficiency  of  current  capabilities  in  modelling  large- 
deformation  failure.  There  presently  does  not  seem  any  obvious  and  rational 
way  of  accounting  for  mechanisms  of  this  type  while  retaining  the 
computational  simplicity  of  the  distributed /imbedded  approach. 


The  lack  of  invariance  of  the  incremental  orthotropic  tangent  moduli 
approach  precludes  its  use  in  a  complex,  dynamic  loading  environment. 
Argument  for  its  use  despite  its  defect  are  simply  not  convincing.  It  would 
seem  that  the  framework  provided  by  anisotropic  damage  mechanics  would 
be  appropriate  for  formulating  invariant  constitutive  equations  capable  of 
representing  cracking.  The  work  of  Resende  [21]  is  a  step  in  this  direction, 
although  it  is  limited  to  isotropic  mechanisms. 


1.  The  cap  model,  in  its  present  form  [15],  does  not  adequately  represent  the 
tension  failure  surface,  tension  stiffening,  anisotropic  elastic  damage  due 
to  cracking,  and  crushing.  Each  of  these  features  can  and  should  be 
improved.  It  should  be  recalled  that  the  cap  model  has  typically  been 
employed  in  highly-confined,  three-dimensional  situations  and  thus  the 
tensile  regime  has  never  been  an  issue.  This  is  not  the  situation  for  the 
plate  ,  slab-like  structures  of  NCEL  and  these  features  need  to  be  addressed 
and  improved. 

2.  Mesh-dependent  softening  moduli  designed  to  desensitize  model 
response  to  degree  of  mesh  refinement  is  an  expedient  methodology,  but 
one  lacking  a  satisfactory  fundamental  basis.  It  is  also  not  clear  how  to  use 
ideas  of  this  kind  when  more  complex  models  are  utilized  (e.g. 
elastoplastic  damage).  An  investigation  into  a  more  fundamental  method 
of  developing  well-posed  models  is  called  for.  In  the  meantime,  mesh- 
dependent  softening  moduli  may  be  employed,  but  it  should  not  be  lost 
sight  of  that  it  is  an  ad  hoc  technology. 

3.  Damage  mechanics  should  be  used  as  a  framework  for  constitutive  theory 
development  so  that  invariant  models  may  be  designed  to  replace  in-place 
noninvariant  models  based  on  incremental  orthotropic  elastic  moduli. 


4.  The  incorporation  of  bond-slip  in  distributed /imbedded  reinforcement 
models  requires  a  breakthrough  that  does  not  seem  imminent.  Current 
developments  should  employ  and  be  content  with  the  perfect  bond 
hypothesis  until  new  ideas  are  forthcoming. 

5.  It  is  probably  worthwhile  to  pursue  a  plane-stress,  plate-shell  model,  as 
well  as  a  more  detailed  three-dimensional  model.  The  former  model  will 
be  more  useful  for  large-scale  structural  modelling  ,  and  the  latter  for 
qualifying  the  former  and  for  detailed  modelling  (e.g.,  in  the  vicinity  of 
supports).  A  new  plane-stress  reinforced  concrete  model  should  be 
developed.  The  existing  cap  should  be  used  as  a  basis  for  further 
development  for  the  three-dimensional  constitutive  model. 

6.  Automated  arc-length  and  time-step  selection  strategies  will  need  to  be 
developed  in  order  to  efficiently  obtain  solutions  to  reinforced-concrete 
structural  problems.  Implicit,  unconditionally  stable  constitutive 
algorithms  should  be  developed  to  circumvent  sub-incremental  stress 
point  stability  limitations. 
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Conclusion 

The  modelling  of  reinforced-concrete  structures  under  severe  dynamic 
loading  represents  a  formidable  challenge.  Many  of  the  key  aspects  of 
reinforced  concrete  are  poorly  understood,  as  evidenced  by  the  numerous 
shortcomings  noted  in  contemporary  models.  Clearly,  one  needs  to  adopt  a 
pragmatic  attitude  and  develop  capabilities  as  best  one  can  in  the  present  and 
immediate  future.  Such  capabilities  will  represent  a  valuable  tool  for  analysis 
and  design  of  critical  structures.  At  the  same  time  one  must  not  be  deceived 
by  the  existence  of  impressive  numerical  capabilities  built  upon  a  foundation 
lacking  in  fundamental  physical  and  mathematical  understanding.  For  the 
foreseeable  future,  reinforced-concrete  modelling  capabilities  will  need  to  be 
exercised  with  extreme  care,  caution  and  insight  to  prevent  naive  and 
erroneous  engineering  judgements  from  being  made.  Pragmatic  and 
fundamental  research  will  need  to  proceed  side  by  side  for  significant  progress 
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